c  
c  solution of Laplace's eq on rectangle
c           0 < x < xmax, 0 < y < ymax
c  
c  computes data for matlab program pclape55
c  error determined using exact solution
c
c	f95 laplace4.f	a.out
c  
c  cgm used to solve linear system
c  
	implicit real*8(a-h,o-z)

	open(unit=7,file="iterations.txt")
	write(7,201)  
  201 format('cgm')
c  201 format(7x,'n',4x,'iterations (cgm)',5x,'time (cgm)')
  
	imin=3
	imax=6
	nmax=15

	ap=0.5*(imax-imin)/(nmax-1.0)
	bp=0.5*(nmax*imin-imax)/(nmax-1.0)
	
	do 80 inn=1,nmax
	pow2=ap*inn+bp
	nx=nint(10**pow2)
	ny=nx
	n=nx*ny
	
	call loop(inn,nx,ny,n,it,xmin_time)
	
c  
c  format statements
c  
	write(7,202)  n, it, xmin_time
  202 format(4x,i8,6x,i4,12x,e13.5)
   80 continue
	close(unit=7)
	end
	
c  
c  main loop 	
c
	subroutine loop(inn,nx,ny,n,it,xmin_time)
	implicit real*8(a-h,o-z) 
	real*4 time1, time2
	parameter(xmax=1.0,ymax=1.0,tol=0.000001)
	dimension b(n),v(n),q(n)
	dimension r(n),d(n),sol(n)

	dx=xmax/(nx+1)
	dy=ymax/(ny+1)
	bb=dy*dy/(dx*dx)
	
	call fourier(n,nx,ny,dx,dy,sol)
	call mult(n,nx,ny,bb,sol,b)
	
	iloops=10
	if(n.gt.1e5)  iloops=5
	
	do 70 in=1,iloops
	
	call cpu_time( time1 )
c
c  define starting value for v
c  
	v=0.0
	call mult(n,nx,ny,bb,v,q)
	r=b-q
	d=r
	rr=dot_product(r,r)
c  
c  cgm iteration
c  
	do 20 it=1,5000
		call mult(n,nx,ny,bb,d,q)
		dad=dot_product(d,q)
		alpha=rr/dAd
		v=v+alpha*d
		r=r-alpha*q
		rr0=rr
		rr=dot_product(r,r)
c		err1=maxval(abs(alpha*d))
		err1=maxval(abs(sol-v))
		if(err1.lt.tol) go to 22
		beta=rr/rr0
		d=r+beta*d
   20   continue
   22	call cpu_time( time2 )
	elapsed_time=time2-time1
	if(in.eq.1)  xmin_time=elapsed_time
	if(in.ne.1.and.elapsed_time.lt.xmin_time) xmin_time=elapsed_time
   70 continue
	return
	end

c  
c  calculate q = Ap
c  
	subroutine mult(n,nx,ny,bb,p,q)
	implicit real*8(a-h,o-z)
	dimension q(n),p(n)
	np=n-nx
	q=2.0*(1.0+bb)*p
	do 10 i=1,n
        	if(i.gt.1) q(i)=q(i)-bb*p(i-1)
        	if(i.lt.n) q(i)=q(i)-bb*p(i+1)
        	if(i.le.np) q(i)=q(i)-p(i+nx)
        	if(i.gt.nx) q(i)=q(i)-p(i-nx)
   10   continue
	do 20 j=1,ny-1
		i=j*nx
		q(i)=q(i)+bb*p(i+1)
		q(i+1)=q(i+1)+bb*p(i)
   20	continue	
	return
	end
c  
c  calculate b
c  
	subroutine right(n,nx,ny,dx,dy,bb,b)
	implicit real*8(a-h,o-z)
	dimension b(n)
	b=0.0
	do 2 ix=1,nx
        	x=ix*dx
        	b(ix)=b(ix)+gB(x)
		ixx=nx*(ny-1)+ix
		b(ixx)=b(ixx)+gT(x)
    2	continue
	do 4 iy=1,ny
        	y=iy*dy
		iyy=nx*(iy-1)
		b(iyy+1)=b(iyy+1)+bb*gL(y)
		b(iyy+nx)=b(iyy+nx)+bb*gR(y)
    4	continue
	return
	end     
c  
c  specify boundary conditions
c      
      function gT(x)
      implicit real*8(a-h,o-z)
      gT=x*(1-x)*(0.8-x)*exp(6*x)
      return
      end
      
      function gB(x)
      implicit real*8(a-h,o-z)
      gB=0
      return
      end
      
      function gR(y)
      implicit real*8(a-h,o-z)
      gR=0.0
      return
      end
      
      function gL(y)
      implicit real*8(a-h,o-z)
      gL=0.0
      return
      end
c  
c  calculate exact solution
c
	subroutine fourier(n,nx,ny,dx,dy,v)
	implicit real*8(a-h,o-z)
	dimension v(n),an(n)
	modes=100
	pi=acos(-1.0)
	do 10 in=1,modes
		a1=cos(0.75*in*pi) - cos(0.25*in*pi)
		an(in)=-2.0*a1/(in*pi*sinh(in*pi))
   10	continue
	do 20 j=1,ny
		y=j*dy
		do 18 i=1,nx
			x=dx*i
			s=0
			do 16 jn=1,modes
				s=s+an(jn)*sinh(jn*pi*y)*sin(jn*pi*x)
   16			continue
			ll=(j-1)*nx+i
			v(ll)=s
   18		continue
   20	continue
	return
	end